clear all;
% clc;

addpath('patient_data');

% ----------------------- Analysis Parameters -----------------------------

pol = -1;
noaxons = 200;

patientgeom = load('femspine_patientgeom.txt');
cordgeom = load('patient_cordgeom.txt');

epidur  = 0; % 0 = intradural, 1 = epidural
subloc  = 1; % 1 = 1 mm above cord; 2 = 1 mm below dura
th_elec = 0; % 0, -10, -20

patients = 1:5;
nopat = length(patients);

dc_th   = zeros(noaxons,nopat);
dr_th   = zeros(noaxons,nopat);
ra_elec = zeros(1,nopat);

% experimental data
isens_exp_sub = [0.20,0.35,0.20,0.55,0.80];
ipain_exp_sub = [0.50,0.50,0.80,1.20,3.10];
isens_exp_epi = [1.30,19.0,7.00,3.00,9.00];
ipain_exp_epi = [2.20,24.5,19.0,6.00,19.0];    

% ----------------------- Importing Data ----------------------------------

model_dcth_epi=cell(nopat,1);
model_drth_epi=cell(nopat,1);
model_iunit_epi=zeros(nopat,1);

model_dcth_sub=cell(nopat,1);
model_drth_sub=cell(nopat,1);
model_iunit_sub=zeros(nopat,1);

diamval=[3,6,9,12,15];
nodiam=length(diamval);

perdcsens_epi=zeros(nodiam,nopat);
perdcsens_sub=zeros(nodiam,nopat);
perdcpain_epi=zeros(nodiam,nopat);
perdcpain_sub=zeros(nodiam,nopat);

for ii = 1:nopat
    
    ptag=['p',num2str(ii),'_'];
    
    model_dcth_epi{ii}=zeros(noaxons,nodiam);
    model_drth_epi{ii}=zeros(noaxons,nodiam);
    model_dcth_sub{ii}=zeros(noaxons,nodiam);
    model_drth_sub{ii}=zeros(noaxons,nodiam);
    
    cd(['p',num2str(ii),'_mrgpotentials']);
    model_iunit_epi(ii)=load([ptag,'epi_bp_current.txt']);
    model_iunit_sub(ii)=load([ptag,'sub_bp_current.txt']);
    cd ..;
    
    cd(['p',num2str(ii),'_swthresholds']);
    for jj=1:nodiam
    
        dcepi_tmp=load([ptag,'epi_bp_dc_act_',num2str(diamval(jj)),'um.txt']);
        
        drepi_tmp=load([ptag,'epi_bp_dc_act_',num2str(diamval(jj)),'um.txt']);
        dcsub_tmp=load([ptag,'sub_bp_dc_act_',num2str(diamval(jj)),'um.txt']);
        drsub_tmp=load([ptag,'sub_bp_dc_act_',num2str(diamval(jj)),'um.txt']);
        
        model_dcth_epi{ii}(:,jj)=abs(dcepi_tmp(:,1))/(2/model_iunit_epi(ii)/1e3);
        model_drth_epi{ii}(:,jj)=abs(drepi_tmp(:,1))/(2/model_iunit_epi(ii)/1e3);
        model_dcth_sub{ii}(:,jj)=abs(dcsub_tmp(:,1))/(2/model_iunit_sub(ii)/1e3);
        model_drth_sub{ii}(:,jj)=abs(drsub_tmp(:,1))/(2/model_iunit_sub(ii)/1e3);
        
        perdcsens_epi(jj,ii)=100*sum(model_dcth_epi{ii}(:,jj)<=isens_exp_epi(ii))/noaxons;
        perdcsens_sub(jj,ii)=100*sum(model_dcth_sub{ii}(:,jj)<=isens_exp_sub(ii))/noaxons;
        perdcpain_epi(jj,ii)=100*sum(model_dcth_epi{ii}(:,jj)<=ipain_exp_epi(ii))/noaxons;
        perdcpain_sub(jj,ii)=100*sum(model_dcth_sub{ii}(:,jj)<=ipain_exp_sub(ii))/noaxons;    
        
    end
    cd ..;
    
end

% dlmwrite('perdc_sens_sub.txt',perdcsens_sub);
% dlmwrite('perdc_sens_epi.txt',perdcsens_epi);
% dlmwrite('perdc_disc_sub.txt',perdcpain_sub);
% dlmwrite('perdc_disc_epi.txt',perdcpain_epi);

% figure;
% subplot(2,1,1),h1=bar(perdcsens_epi');
% set(gca,'XTick',[],'YTick',25:25:100,'FontSize',30);
% title('Extradural','FontSize',30);
% xlim([0.5,5.5]);
% ylim([0,105]);
% 
% subplot(2,1,2),h2=bar(perdcsens_sub');
% set(gca,'YTick',25:25:100,'FontSize',30);
% title('Intradural','FontSize',30);
% xlim([0.5,5.5]);
% ylim([0,105]);
% set(h1,'LineWidth',3);
% set(h2,'LineWidth',3);
% colormap('gray');
% legend('3 \mum','6 \mum','9 \mum','12 \mum','15 \mum');
% 
% figure;
% subplot(2,1,1),h3=bar(perdcpain_epi');
% set(gca,'XTick',[],'YTick',25:25:100,'FontSize',30);
% title('Extradural','FontSize',30);
% xlim([0.5,5.5]);
% ylim([0,105]);
% 
% subplot(2,1,2),h4=bar(perdcpain_sub');
% set(gca,'YTick',25:25:100,'FontSize',30);
% title('Intradural','FontSize',30);
% xlim([0.5,5.5]);
% ylim([0,105]);
% set(h3,'LineWidth',3);
% set(h4,'LineWidth',3);
% colormap('gray');
% legend('3 \mum','6 \mum','9 \mum','12 \mum','15 \mum');

% ----------------------- Model Thresholds --------------------------------

isens_model_epiLL=zeros(nopat,1);
isens_model_epiLU=zeros(nopat,1);
isens_model_epiUL=zeros(nopat,1);
isens_model_epiUU=zeros(nopat,1);

isens_model_subLL=zeros(nopat,1);
isens_model_subLU=zeros(nopat,1);
isens_model_subUL=zeros(nopat,1);
isens_model_subUU=zeros(nopat,1);

ipain_model_epiLL=zeros(nopat,1);
ipain_model_epiLU=zeros(nopat,1);
ipain_model_epiUL=zeros(nopat,1);
ipain_model_epiUU=zeros(nopat,1);

ipain_model_subLL=zeros(nopat,1);
ipain_model_subLU=zeros(nopat,1);
ipain_model_subUL=zeros(nopat,1);
ipain_model_subUU=zeros(nopat,1);

for ii=1:nopat
    
    isens_model_epiLL(ii)=min(min(model_dcth_epi{ii}));
    isens_model_epiLU(ii)=max(min(model_dcth_epi{ii}));
    isens_model_epiUL(ii)=min(max(model_dcth_epi{ii}));
    isens_model_epiUU(ii)=max(max(model_dcth_epi{ii}));
    
    isens_model_subLL(ii)=min(min(model_dcth_sub{ii}));
    isens_model_subLU(ii)=max(min(model_dcth_sub{ii})); 
    isens_model_subUL(ii)=min(max(model_dcth_sub{ii}));
    isens_model_subUU(ii)=max(max(model_dcth_sub{ii})); 
    
    ipain_model_epiLL(ii)=min(min(model_drth_epi{ii}));
    ipain_model_epiLU(ii)=max(min(model_drth_epi{ii}));
    ipain_model_epiUL(ii)=min(max(model_drth_epi{ii}));
    ipain_model_epiUU(ii)=max(max(model_drth_epi{ii}));
    
    ipain_model_subLL(ii)=min(min(model_drth_sub{ii}));
    ipain_model_subLU(ii)=max(min(model_drth_sub{ii})); 
    ipain_model_subUL(ii)=min(max(model_drth_sub{ii}));
    ipain_model_subUU(ii)=max(max(model_drth_sub{ii})); 
    
end

% ----------------------- Processing Data ---------------------------------

figure;

% EPIDURAL
x=1:nopat;
sepiL=isens_model_epiLL;
sepiU=isens_model_epiUU;
pepiL=ipain_model_epiLL;
pepiU=ipain_model_epiUU;

subplot(1,2,1);
hold on;
plot(patients-0.2,isens_exp_epi,'ko','LineWidth',4,'MarkerSize',15);
for kk=1:nopat
    plot([x(kk)-0.05,x(kk)-0.05],[sepiL(kk),sepiU(kk)],'k-','LineWidth',4);
    plot([x(kk)+0.05,x(kk)+0.05],[pepiL(kk),pepiU(kk)],'k--','LineWidth',4);
end
plot(patients+0.2,ipain_exp_epi,'ko','LineWidth',4,'MarkerSize',15,...
     'MarkerFaceColor','k');

hold off;
set(gca,'FontSize',36,'XTick',1:5,'YScale','Log','YTick',[0.01,0.1,1,10,50]);
xlim([0.5,5.5]);
ylim([0.1,75]);
title('Extradural','FontSize',36);
xlabel('Patient','FontSize',36);
ylabel('Stimulation Amplitude (mA)','FontSize',36);
axis square;

% SUBDURAL
ssubL=isens_model_subLL;
ssubU=isens_model_subUU;
psubL=ipain_model_subLL;
psubU=ipain_model_subUU;

subplot(1,2,2);
hold on;
plot(patients-0.2,isens_exp_sub,'ko','LineWidth',4,'MarkerSize',15);
for kk=1:nopat
    plot([x(kk)-0.05,x(kk)-0.05],[ssubL(kk),ssubU(kk)],'k-','LineWidth',4);
    plot([x(kk)+0.05,x(kk)+0.05],[psubL(kk),psubU(kk)],'k--','LineWidth',4);
end
plot(patients+0.2,ipain_exp_sub,'ko','LineWidth',4,'MarkerSize',15,...
     'MarkerFaceColor','k');
set(gca,'FontSize',36,'XTick',1:5,'YScale','Log','YTick',[]);
hold off;
set(gca,'FontSize',36,'XTick',1:5,'YScale','Log');
title('Intradural','FontSize',36);
xlabel('Patient','FontSize',36);
xlim([0.5,5.5]);
ylim([0.1,75]);
axis square;

